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We derive an efficient method for the insertion of structured particles in grand canonical Monte 
Carlo simulations of adsorption in very confining geometries. We extend this method to path integral 
simulations and use it to calculate the isotherm of adsorption of hydrogen isotopes in narrow carbon 
nanotubes (2D confinement) and slit pores (ID confinement) at the temperatures of 20 K and 77 K, 
discussing its efficiency by comparison to the standard path integral grand canonical Monte Carlo 
algorithm. 

We use this algorithm to perform multicomponent simulations in order to calculate the hydrogen 
isotope selectivity for adsorption in narrow carbon nanotubes and slit pores at finite pressures. 

The algorithm described here can be applied to the study of adsorption of real oligomers and 
polymers in narrow pores and channels. 



I. INTRODUCTION 

The grand canonical (GC) Monte Carlo (MC) method 
is an efficient way to calculate thermodynamic proper- 
ties of adsorbed fluids via computer simulations It 
has recently been extended to take into account quan- 
tum diffraction effects using the path integral (PI) for- 
mulation of quantum mechanics if the adsorption of light 
particles, such as helium atoms or hydrogen molecules, is 
investigated at low temperatures^i^ 

Use of the method for the study of hydrogen adsorp- 
tion in very confining geometries has led to the discovery 
of a very strong isotope effect 1^1^ The adsorption of hy- 
drogen isotopes in carbon nanotubes has been thoroughly 
investigated as a function of temperature and nanotube 
radius, and it has been found that these nanostructured 
materials might be used, in principle, to obtain a very 
efficient isotope separation, a phenomenon referred to as 
"quantum sieving" and studied theoretically for a long 
timcj^i^iiiSiiSii^'iiii^'i^ii^ Rcccntly the first experimental 
confirmations have been reportedJ^ii^ 

From a simulation point of view the PIMC method 
maps the statistical mechanics of N quantum particles 
into the statistical mechanics of N ring polymers, each 
formed by P particles. The mapping is exact in the limit 
P ^ (XI, even though, according to the temperature in- 
vestigated, convergence of the results can be obtained 
when P is of the order of 10-100. 

The standard path integral grand canonical 
(PIGCMC) simulation procedure was developed by 
Johnson and collaborators^ who used ring polymer 
configurations drawn from an ideal gas distribution as 
insertion candidates. While this method is effective for 
studying adsorption in large pores, it is very inefficient 
when applied to very narrow geometries. The reason of 
this inefficiency stems from the fact that in a strongly 
confined system at low temperatures the average kinetic 
energy of the adsorbed particles might be much higher 
than the average thermal kinetic energy of a free particle. 
In the path integral formulation this quantity is related 
to the size of the ring polymers: polymers are smaller - 
on average - when they represent particles with higher 



kinetic energy. As a consequence, if one tries to insert 
ring polymers drawn from a ideal gas distribution at the 
desired temperature in a strongly confining geometry, 
the likelihood of a rejection is very high, because the 
"average" ideal gas configuration just does not fit in the 
available space within the adsorbent. 

This phenomenon has been recognized in the early 
PIGCMC simulations, since very low acceptance ratios 
were observed when trying to insert ideal gas ring poly- 
mers in very narrow carbon nanotubesJ^ The same phe- 
nomenon has also been observed by independent studies 
in the case of adsorption in carbon slit pores and and has 
been called "quantum polymer shrinking" by Kowalczyk 
et alA and it is a consequence of the well known fact 
that more confining geometries enhance the zero point 
energy of adsorbed particles. This effect is even more 
pronounced when the rotational degrees of freedom of 
the hydrogen molecule are taken into account, leading to 
a dramatic enhancement of the sieving propertiesj^^^iii 

We propose a new algorithm for insertion moves in 
classical and quantum (path integral) grand canonical 
Monte Carlo simulations, which is particularly suitable 
in the case of strong confinement. We show that the ac- 
ceptance ratio of insertion moves can be up to two orders 
of magnitude higher than the one obtained using ideal gas 
ring polymers as trial configurations. As a consequence 
the equilibration and production runs of a PIGCMC sim- 
ulation can be made correspondingly shorter, consider- 
ably reducing the time needed for the calculations. 

This method, which we call Boltzmann bias grand 
canonical Monte Carlo, uses as trial configurations ring 
polymers obtained by an independent path integral 
Monte Carlo simulation of particles which are adsorbed 
within the same adsorbent, but do not interact between 
themselves, so that the "shape" of the ring polymers to 
be inserted is already the correct one and consequently 
the insertion candidates are more likely to fit within the 
confining geometry. 

In the following we derive the expressions for the trial 
and acceptance probability, both in the classical and in 
the quantum case. Although the Boltzmann bias method 
is not likely to produce any improvement in classical sim- 
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ulations of structureless particles with respect to other 
existing algorithms, we present the derivation in detail 
for the classical case. This will be done to fix the notation 
before passing on to the quantum case. We also present 
this algorithm in the case of adsorption of oligomers and 
polymers within confining geometries. 

We will apply the Boltzmann bias method to the cal- 
culation of adsorption isotherms and finite pressure se- 
Icctivities for mixtures of hydrogen isotopes in narrow 
carbon nanotubes and slit pores, and compare our re- 
sults with those obtained using the standard PIGCMC 
method. 

The analysis of the efficiency of the Boltzmann bias 
method will show that one can obtain acceptance ratios 
up to four orders of magnitude higher than the ones ob- 
tained using the standard insertion algorithm, dramati- 
cally reducing the amount of time needed to perform the 
computer simulations. 



II. BOLTZMANN BIAS METHOD 



Classical simulations 



moves that change the number of particles in the sys- 
tem. It is usually sufficient to consider moves in which 
the number of particles is changed by unity. 

The expression for the transition probability Wi^j of 
the Markov chain between states i and j having and 
AT -|- 1 particles respectively can be obtained by applying, 
as usual, the detailed balance condition, i.e. 



Wn- 



Wn+1 



e-^^^d'x^+, (2) 



where AU = Un+i — Un is the difference between the 

potential energy of the two states i and j. The transi- 
tion probability W is written as the product of a trial 
probability T and an acceptance probability A. In the 
standard grand canonical Monte Carlo algorithm the trial 
probability for an insertion is such that a test particle is 
inserted uniformly within the volume and the the prob- 
ability of a removal is just one, hence 



d^XN 



+1 



V 



(3) 



The GCMC method is based on the generation of a 
Markov chain having as equilibrium probability density 
the quantity 



p(<'i)(a;i,...,arAr) = i 



(1) 

where A'' is the variable number of particles, and V, T 
and /z are the system's volume, the temperature and 
the chemical potential, respectively, which are all kept 
fixed. The quantity A = h/y^2TrmkBT is the ther- 
mal de Broglie wavelength of the particles of mass m 
and S 'J2n ^^^^Qi^i^T^^) grand canonical par- 

tition function. Qn{V,T) is the canonical partition 
function, /3 — l/ksT and ks is the Boltzmann constant. 
Finally, U{xi, . . . ,Xn) is the interaction potential energy 
between the N particles and/or between the particles and 
an external potential (which can model, for example, an 
adsorbent). 

The Markov chain is built by moves that either keep 

the number of particles constant (and for which the stan- 
dard canonical Monte Carlo algorithms can be used) and 



so that the ratio between the acceptance probabilities 
corresponding to the insertion into and removal from the 
system of a particle is 



V 



An+i^n {N + l)Af 



(4) 



a condition that can be satisfied, e.g., by the standard 
Metropolis choice 



= min 

An+i^n = min 



V 



' {N + l)Al 
1, ^ T-r^ — e '"''e 



(7V + 1)A3 ' 



V 



(5) 



The Boltzmann bias methods stems now from the fact 
that in many cases the interaction potential energy can 
be written as the sum of the potential energy of interac- 
tion with the fluid particles between themselves (Uff) 
and the interaction potential energy with a substrate 
(Usf), the latter being additive in the number of par- 
ticles, i.e. 



N 



U{X1, . . .,Xn) = Uff{xi, . . .,Xn) + Usf{xi, . . .,Xn) = Uff{xi,. . .,Xn) + '^USF{Xk) (6) 
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so that one can devise an insertion move such that the solid-fluid interaction potential of the particle to be in- 
probability of insertion within the volume V is not uni- serted, exp [— /3usF(a;jv-i-i)]- One can then write the trial 
form, but proportional to the Boltzmann factor of the 
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probability for an insertion as 



exp [-/3MSF(a:Ar+i)] cflxN+1 _ exp[-PusF{xN+l)] 



exp [-^UsF(a;Ar+l)] d^XN+i 



(7) 



r 



where the last equality defines the quantity p.. With this 
choice for the trial probability of an insertion (leaving 
Tn+i^n = 1), one gets for the acceptance probability 
the ratio 



V 



An+i^n {N + 1)A3 



g/3(A'-A')g-A(7FF(a;i, 



(8) 



that depends only on the difference of the fluid-fluid po- 
tential energy between the original configuration and the 
trial one. 

One way to insert particles with the trial probability 
given by Eq. ([7]) is to perforin, in parallel with the grand 
canonical simulation, a standard canonical Monte Carlo 
simulation of some particles within a similar box and 
at the same temperature. In this second simulation the 
particles interact only with the substrate and not among 
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N+l 



themselves, so that at equilibrium they sample the distri- 
bution given by Eq. ([7]). One can then get the coordinates 
for the GCMC trial insertions by drawing particles' po- 
sitions from the configurations generated by this second 
simulation. 

In the case of classical adsorption of structureless par- 
ticles this method is not particularly more efficient than 
other standard methods based on the insertion in the 
subvolume of the simulation cell where the solid-fluid 
potential energy is less than a given threshold. This 
free volume Vfreo can be easily pre-calculated and used 
in the standard grand canonical algorithm for insertions, 
Eq. ([3]), instead of the whole cell volume V. 

In the case of classical adsorption of molecules, the 
detailed balance condition ^ can be written as 



(9) 



where x'^i denotes the center of mass position and ^ 
are the variables describing the intra-molecular degrees 
of freedom (such as Euler angles in the case of small rigid 
molecules or torsion angles for small oligomers). The 
quantity ?7int(0 is just the internal energy corresponding 
to the configuration given by the coordinates ^. One can 
see from Eq. ([9]) that the transition probability depends 
on the internal state of the molecule ^ as well as the po- 
sition of the center of mass for the trial insertion. In 
the limit of low density the fluid-fluid contribution to the 
acceptance probability can be neglected and this quan- 
tity is mainly determined by the solid-fluid part of the 
interaction. In the case of adsorption in narrow geome- 
tries, many of the trial configurations f for a given value 
of the center of mass position x'^^ are likely to be re- 
jected, since they would result in strong overlaps of parts 
of the molecule with the adsorbent. In this case a much 
greater acceptance ratio could be obtained if the config- 
urations ^ were chosen among the ones most likely to fit 
at a given position x'^-^^. The Boltzmann bias method 
proposed here is based on choosing the the trial configu- 
rations with a probability proportional to the Boltzmann 



factor of the molecule to be inserted at the position x'^-^, 
i.e. 

By writing the normalization factor as 

Vel^P = J e-/3^sp(.5^-^,.«)e-/3C/.„.(«) d^xf^^d^ (11) 

one obtains for the acceptance probability an expression 
analogous to the one reported in Eq. ([5]). The trial con- 
figurations can be easily generated by running an NVT 
simulation of molecules within the adsorbent, with the 
fiuid-fiuid interaction switched off. 

Since the solid-fiuid interaction has already been taken 
into consideration by the choice of the trial probability, 
the acceptance probability now depends on the values of 
the fiuid-fluid interaction between the molecule to be in- 
serted and the ones already present within the simulation 
box, as well as on the value of the renormalized chemical 
potential fJ. — p.- One can then expect a very high accep- 
tance ratio for all the thermodynamic state points where 
the density of the adsorbed molecules is not liquid-like. 
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B. Path integral simulations 

The path integral grand canonical method is based on 
the Trotter factorization of the many-body Hamiltonian. 
The quantum partition function of N interacting parti- 
cles can be written as 



1 



P.) 



(12) 



where T and V are the kinetic and potential energy op- 
erators, respectively, and the expectation value is taken 
over a complete set of N particle states. The sum is over 
all the permutations tt of the particles, represented in 
the Hilbert space by the operator Pn- which also takes 
into account the bosonic or fermionic nature of the par- 
ticles. In the following we will not consider effects due 
to quantum statistics and hence we will approximate the 
sum with its leading term, corresponding to the identity 
permutation. This is a good approximation as long as 



the thermal de Broglie wavelength A does not exceed the 
hard core radius of the particles under consideration.— 
By using the Trotter formula 



= lim 



(13) 



with a finite but large enough Trotter index P, the quan- 
tum partition function of N particles can be mapped into 
the classical partition function of NP particles, such that 
each original particle corresponds to a ring polymer of P 
"beads" The beads within each polymer interact with 
their two nearest neighbors via an harmonic interaction, 
and we will call this part of the interaction "internal". 
Denoting by Xi = {x^^ , . . . , a;-^-*} the 3P coordinates of 
the beads of the ring polymer corresponding to the i-th 
particle, and applying the same reasoning as described 
before in the classical case, one obtains for the path in- 
tegral case 



J 



1 / p3P/2 



Wn+i^n N+1\ A^p 



„/3Aig-/3A(7FF(Xi,...,Xjv+i)g-/3(7sF(XN+i)g-/3;7i„t(Xjv + i) d^Xw+l 



(14) 



r 



where we have defined 

p N 

U^{X,,...,X^) = pJ2T. "FF(N?-xf 1015) 



p—1 i<j—l 

p 



Usf{Xn+i) = ^^^^SFix^N+l) 



Uint{XN+l) = y ^ pXr+i - a::]v+i 



JP+1) 



(16) 



(17) 



K 



27rP/A2, 

, Xn^ 



i^UFFiXi, . . . ,Xn+i) — 

Uff{Xi, . . . , Xn) and we have 



Uff{Xi, . . . , ^N+l) 

used the convention x']^^^'' = xj^l^j^. In the path 
integral GCMC case the ratio between the transition 
probabilities connecting states with different number of 
particles depends also on the internal state of the added 
polymer and it is customary to draw internal states from 
an ideal gas path integral simulation.^ 

This way of proceeding is highly inefficient in the case 
of strong confinement, where one expects that the distri- 
bution of shapes for the adsorbed ring polymers, having 
to account for a high kinetic energy due to Hciscnberg 
indetermination, will be considerably different from the 
equilibrium distribution of an ideal gas. In particular it 



turns out that polymers corresponding to states of higher 
kinetic energy tend to have a smaller size.^' As a conse- 
quence the probability of rejecting an insertion move is 
particularly high, because almost all of the ideal gas ring 
polymer configurations do not easily fit into the available 
space within the adsorbent. 

In this case the Boltzmann bias method is likely to im- 
prove the situation. Instead of drawing a ring polymer 
configuration from an ideal gas distribution, the trial in- 
sertion probability is chosen to be proportional to the 
probability of finding an isolated ring polymer within the 
adsorbent, i.e. 



/ p3P/2 \ 



By using the change of variables 



(18) 



r 
Ari 



(P-1) 



(19) 



the normalization factor can be written as the integral 
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1 

A3 



-l3UsF{r,Ari,...,Arr-l) 



i^ring(Ari, 



P-1 

.Arp_i) Jl At- 
p=i 



d^r = Ve 



(20) 



where Fi.ing(Ari, . . . , Arp_i) is the probabiUty of having 
a configuration of an ideal gas ring polymer with sep- 
arations Arp between the neighboring beads, and it is 
derived in Appendix El The last equality defines the 
quantity /i analogously to the classical definition given in 
Eq. Q. 

With these definitions the ratio of the acceptance prob- 
abilities for moves which change the number of particles 
in the path integral case is 



iAr->iV+i 



V 



g/9(M-p)g-A;7FF(Xi,...,Xjv + i) 



(21) 

completely analogous to the classical case of Eq. By 
"undoing" the Trotter factorization in Eq. ([20)1 one can 
show that 



Ve 



A' 



exp 



P 

2m 



-fJEi 



(22) 

where Ei are the quantum energy levels of a single par- 
ticle interacting with the adsorbent. 

The Boltzmann bias PIGCMC method is particularly 
efhcient in the case where adsorption is simulated in ma- 
terials having an (approximate) continuous translational 
symmetry as it happens, e.g., in carbonaceous materials 
such as carbon nanotubes and slit pores, where it has 
been shown that the corrugation of the potential energy 
can be neglected without sacrificing accuracy.^ In this 
case the ring polymer configurations for trial insertions 
can be further translated by a random value along the 
direction(s) of translational symmetry. This prescription 
avoids the need to wait for the polymers in the zero- 
pressure simulation box to diffuse away from positions 
already used for trial insertions. 

In the case of simulations in geometries that do not 
possess continuous translational symmetry, as it might 
happen if the effect of nanotube corrugation on adsorp- 
tion had to be addressed or in the case of adsorption in 
amorphous materials, one should perform long stages of 
NVT MC in the zero-pressure box until a new configu- 
ration, uncorrelated with the previous one, is generated. 
This is necessary in order for the new polymers used for 
trial insertions not to end up on top of polymers already 
present in the simulation box, resulting in a certain re- 
jection. 

In the general case another method for generating trial 
polymers, which seems particularly suited for path inte- 
gral simulations, could be used. One starts from con- 
figurations where N independent classical point particles 
are adsorbed at the target temperature T (which can 



be generated very efficiently, even in a parallel environ- 
ment, either by GCMC or molecular dynamics), and use 
them to build N ring polymers where all the beads are 
at the position of the corresponding classical particles. 
These polymers can be efficiently equilibrated using hy- 
brid MC^'^" and then used as insertion candidates. This 
approach might be more efficient than the direct one out- 
lined above, and its performance will be analyzed in fu- 
ture studies. 

There is another issue to be addressed in the case of 
simulation of adsorption in more general geometries, i.e. 
the calculation of the parameter /2. In a general geometry 
the calculation of the energy levels of a single particle, 
needed to compute /I according to Eq. (|22[) . can be very 
computationally demanding and in this case it could be 
worth it using directly Eq. ([^0]) . In fact the left hand side 
of Eq. (|20p consists of a three dimensional integration 
over the variable r of the e"'^'^^''''"''^'"^' ' ''^'"^-'^'' averaged 
over ring polymer configurations drawn from an ideal gas 
distribution, which can be generated rather efficiently;^ 



III. ADSORPTION IN NARROW CARBON 
NANOTUBES 

A. Pure fluid isotherms 

In order to validate the Boltzmann bias method we 
compare the results with standard path integral grand 
canonical Monte Carlo simulations of adsorption in car- 
bon nanotubes, where insertion moves are performed by 
taking ring polymer configurations from an ideal gas sim- 
ulation performed in parallel with the main computation. 

We model the carbon nanotubes as smooth tubes of 
given geometrical radius R. The solid-fluid interaction 
depends only on the distance of the particle from the tube 
axis, and we used the expression for the average potential 
reported in Refs. [2^l23ll2^ . We chose two values of the 
geometrical radius: R = 3.6 A that corresponds to the 
(2,8) tube and i? = 3.1 A that corresponds to the (3,6) 
tube. 

Earlier path integral simulations of hydrogen adsorp- 
tion in these tubes showed that at the temperature of 
T = 20 K, quantum effects are very important to de- 
scribe the actual behavior of hydrogen adsorbed in the 
(3,6), and moderately important in the case of the (2,8) 
tnhe^^ At a higher temperature of T = 77 K, quantum 
effects on adsorption in the (2,8) tube are expected to be 
smaller. 

In our simulations the fluid-fluid interaction between 
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hydrogen atoms is assumed to be of the Lennard- Jones 
type, with the Buch parameters'^ eHH/^s — 34.2 K and 
cthh = 2.96 A, and the soHd-fluid hydrogen-carbon inter- 
action is calculated using the Lorentz-Berthelot mixing 
rules^ together with the Steele parameters-- for the car- 
bon atom, ecc/ks = 28.0 K and ctcc = 3.4 A. The quan- 
tity /2 appearing in Eq. (|2ip has been calculated using 
Eq. ([^^ . The energy levels of a single hydrogen particle 
interacting with a carbon nanotube have been obtained 
by direct diagonalization of the quantum Hamiltonian, 
using as a basis set the eigenf unctions of a free parti- 
cle confined within a rigid cylinder of radius R. In the 
case of the (3,6) nanotube, the first 183 energy states 
of the free particle were sufficient to reach convergence, 
whereas in the case of the (2,8) tube we obtained conver- 
gence of the results using 226 states. The values of /2 as 
well as other relevant single particle properties of hydro- 
gen isotopes confined in carbon nanotubes are reported 
in Table [J The results of these single particle calcula- 
tions were used to fix the value of the Trotter number P 
as a function of the temperature: we performed path in- 
tegral canonical simulations of hydrogen adsorbed within 
the tubes switching off the interparticle interaction, pro- 
gressively increasing the value of P until we obtained the 
same average values for the kinetic and potential energies 
as those reported in Table [H We reached convergence us- 
ing P = 64 for Ha at T = 20 K and using P = 16 for 
H2 at T = 77 K. The same values of the Trotter index P 
obtained for H2 were used for the other isotopes. 

We show in Fig. [1] the adsorption isotherms obtained 
using both the Boltzmann bias method and the stan- 
dard procedure for path integral grand canonical Monte 
Carlo. All of the simulations have been performed us- 
ing 1.5 X 10^ MC moves per thermodynamic point, after 
7.5 X 10^ moves for equilibration. We consider adsorption 
in an isolated nanotube of height Z = 400 A. Out of 100 
Monte Carlo moves, 80 were insertion/deletion attempts 
and 20 were hybrid MC moves 1^1^ 

A first test on the efficiency of the Boltzmann bias 
method is performed by calculating the adsorption of hy- 
drogen in the (2,8) tube at T = 77 K. The results of the 
Boltzmann bias method are compared with the standard 
prescription for insertion moves in Fig. [Ija), where it is 
shown that the two methods lead to the same results as 
far as adsorption is concerned. 

The results for the same tube at the lowest tempera- 
ture, shown in Fig. [TJb) show a slight disagreement, es- 
pecially in the region around saturation. The isotherms 
obtained with the two methods for the (3,6) tube at 
T = 20 K, reported in Fig. [T]^c), are significantly dif- 
ferent. The origin of the difference is due to the poor 
performance of the standard insertion method, that re- 
sults in a very high rejection rate and hence in a poor 
convergence towards the equilibrium density. The result 
of the standard PIGCMC simulation of H2 adsorption in 
the (3,6) tube using twice as many steps for equilibration 
and production, also reported in Fig. [IJc), confirms this 
observation. As is apparent this second isotherm tends 




Fugacity [bar] 



120| 




Fugacity [bar] 



FIG. 1: Isotherms of adsorption for hydrogen in carbon nan- 
otubes. (a) (2,8) tube at T = 77 K, (b) (2,8) tube at T = 20 K 
and (c) (3,6) tube at T = 20 K. The circles denote the results 
using the Boltzmann bias method, whereas the triangles are 
the results of simulations using ideal gas ring polymers as in- 
sertion candidates. The diamonds in panel (c) are the result 
of a standard PIGCMC calculation with twice as many MC 
steps as the other one. 



to converge to the correct result, but still shows evident 
signs of poor equilibration. 

In order to prove this fact, we show in Fig. [2] the inser- 
tion and deletion acceptance ratios for all the simulations 
we have performed. 

We notice that the insertion and deletion acceptance 
ratios have the same values for each state point, except 
than in the case of the (3,6) tube simulated using the 
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Nanotube type 


Adsorbate 


Temperature (K) 


U/kB (K) 


Kinetic Energy/fcs (K) 


Potential Energy/fe_g (K) 


(2,8) 


Ha 


20 


-1245 


124.0 


-1392 


(2,8) 


Ta 


20 


-1314 


69.1 


-1428 


(2,8) 


Ha 


77 


-1054 


167.4 


-1383 


(2,8) 


Ta 


77 


-1088 


132.5 


-1404 


(3,6) 


Ha 


20 


-281.3 


329.9 


-629.7 


(3,6) 


Ta 


20 


-524.3 


187.2 


-751.0 


(3,6) 


Ha 


77 


-99.05 


358.6 


-628.6 


(3,6) 


Ta 


77 


-281.0 


219.5 


-747.7 



TABLE I: Single particle properties of hydrogen isotopes confined in carbon nanotubes, obtained by direct diagonalization of 
the single particle Hamiltonian. The quantity is defined in Eq. (|22|) . 



standard practice of inserting ideal gas ring polymers. A 
simulation run that achieves the same acceptance ratios 
for insertion and deletion is well equilibrated. The fact 
that this situation is not verified in the (3,6) tube using 
standard moves means that equilibration has not been 
achieved even after 7.5 x 10^ steps. 

Under the same conditions the Boltzmann bias method 
does not show any sign that convergence has not been 
reached, and we have checked that one obtains the same 
isotherm using a larger number of Monte Carlo moves 
(2.5 X 10® for production and half as many for equilibra- 
tion, result not shown). 

We notice that at each thermodynamic point that has 
been investigated the acceptance ratios for the Boltz- 
mann bias method are at least one to two order of mag- 
nitude higher than the ones obtained using the standard 
method, and they can be over four orders of magnitude 
higher than the standard method in very confining ge- 
ometries at low pressure, such as it happens in the case 
of hydrogen adsorption at T = 20 K in (3,6) tube (see 
Fig. He)). 

As might have been expected, the acceptance ratio 
decreases as the density of adsorbed molecules is in- 
creased. This is due to the higher probability of over- 
lap between the particle to be inserted and the particles 
already present in the simulation cell as the density of 
adsorbed molecules is increased. In order to raise the ac- 
ceptance ratio in this case, other kinds of biases can be 
considered, such as the configuration-bias Monte Carlo. 
In any case the acceptance ratios for the Boltzmann bias 
method are one to two orders of magnitude higher than 
the acceptance ratios obtained using ideal gas ring poly- 
mers as candidates. This means that simulations using 
the Boltzmann bias method need one to two order of 
magnitude less steps to achieve equilibration as the stan- 
dard method. 



Nanotube type 


Mixture 


Temperature (K) 


Zero pressure selectivity 


(3,6) 


Ta/Ha 


20 


181000 


(3,6) 


Ta/Ha 


77 


10.5 


(3,6) 


Da/Ha 


20 


5000 


(3,6) 


Da/Ha 


77 


5.5 


(2,8) 


Ta/Ha 


20 


32 


(2,8) 


Ta/Ha 


77 


1.55 


(2,8) 


Da/Ha 


20 


12.5 


(2,8) 


Da/Ha 


77 


1.40 



TABLE II: Zero pressure selectivities of hydrogen mixtures in 
carbon nanotubes 



at low temperatures and under conditions of strong con- 
finement^>^i>^^^>i^^^ 

The differential adsorption of two species A and B is 
measured by the selectivity S{A/B), defined as 



S{A/B) 



Va/vb 



(23) 



where Xi is the molar fraction of specie i in the adsorbed 
phase, and yi is the molar fraction of the same specie in 
the gas phase, under conditions of thermodynamic equi- 
librium. When the pressure is so low that the fluid-fluid 
interaction between the particles can be neglected, the 
selectivity can be written as a function of the single par- 
ticle energy levels of species A and B in the adsorbed 
phasoiS 



( aiE. 



3/2 



(24) 



where rrii are the masses of the two species and Qi is the 
single particle partition function 



B. Zero and finite pressure selectivity in narrow 
tubes 

As is well known quantum diffraction effects result in 
different adsorption properties of various isotopes of the 
same specie, a phenomenon which is particularly evident 



= ^exp(-/34' 



(25) 



It can be seen that, in the limit T ^ 0, the selectivity 
is a function of the difference of the zero point energies 
of the two adsorbed species. Values of the zero pressure 
selectivity for mixtures of hydrogen isotopes adsorbed in 
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FIG. 2: Insertion and deletion acceptance ratios for the pro- 
duction runs of the grand canonical simulations, (a) (2,8) 
tube at T = 77 K, (b) (2,8) tube at T = 20 K and (c) (3,6) 
tube at r = 20 K. The circles denote the acceptance ratios 
obtained using the Boltzmann bias method, whereas the tri- 
angles are the acceptance ratios from simulations using ideal 
gas ring polymers as insertion candidates. The solid lines 
join the points corresponding to the deletion acceptance ra- 
tios (not highlighted by any symbol). 



carbon nanotubes are reported in Tab. [TTl As might be 
expected, the selectivity is higher when the mass ratio be- 
tween the isotopes is higher (all other things being equal) , 
and has a strong temperature dependence. The values re- 
ported there can be compared with the values reported 
in literature for molecular sieves (like zeolites) used in 
gas separation, where typical values are in the range 10- 



FIG. 3: Selectivity for a T2/H2 mixture in the (3,6) carbon 
nanotube at T = 20 K as a function of pressure. The simu- 
lation has been performed by assuming a mole fraction of T2 
in the bulk phase equal to 3/T2 = 5 x 10~^. The horizontal 
dashed line corresponds to the zero pressure value, reported 
in Tab. HIl 

100, depending on the components of the mixture to be 
separatedi22i2£ 

Challa et al. have investigated the effect of finite pres- 
sure on the selectivity, by performing multicomponent 
grand canonical path integral Monte Carlo simulations 
and calculating the selectivity directly from the defini- 
tion of Eq. They noted a progressive inefficiency 
in the acceptance ratios for insertion especially at high 
pressures, and hence were unable to investigate the finite 
pressure selectivity in conditions of strong confinement. 

The use of the Boltzmann biased algorithm allows one 
to gain a factor of 100 in the insertion efficiency, and cor- 
respondingly reduce the time needed for equilibration. 
The results of our simulations for the pressure depen- 
dence of the T2/H2 selectivity in the (3,6) nanotube are 
reported in Fig. [3l We have used nanotubes of different 
length for different pressures: the nanotube lengths were 
adjusted to have between 150 and 300 particles in the sys- 
tem after equilibration. In order to have an almost equal 
amount of hydrogen and tritium within the simulation 
cell we have chosen to perform the simulation by impos- 
ing a bulk mole fraction of T2 equal to — b x 10~^. 
We have used 2 x 10^ steps for equilibration and 4 x lO'^ 
steps for production at each state point. The probability 
of performing an insertion/deletion has been set to 0.8. 

One can see that, according to what has been observed 
in analogous confining geometries, the selectivity actually 
increases from its zero pressure value and tends to satu- 
rate to a value which is some 40% higher that the zero 
pressure one. This behavior can be understood using 
a simple mean field model by assuming that, upon sat- 
uration, the particles are uniformly arranged within the 
tube, so that each particle is confined also in the direction 
of the tube axis. Using the simple theory of Eq. p4)) . the 
selectivity upon saturation can be approximately writ- 
ten as the zero pressure selectivity times the selectivity 
of a one dimensional harmonic oscillator corresponding 
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to the motion along the tube axis. A simple calculation 
shows that one does indeed get the right magnitude of 
the increase.— 



IV. ADSORPTION AND SELECTIVITY IN SLIT 
PORES 

We have also applied our algorithm to the calculation 
of adsorption isotherms and quantum sieving effects of 
hydrogen in a narrow slit pore. We have modeled a 
slit pore as being formed by two graphite planes sepa- 
rated by a distance H along the direction z, using the 
same parameters for carbon-hydrogen interaction as in 
the case of nanotubes described above. We have obtained 
a smooth solid-fluid interaction potential Vsntiz) by av- 
eraging the values of the potential obtained by atomistic 
calculation over planes at constant height z. As a re- 
sult of this approximation the quantum states of a parti- 
cle adsorbed within the slit pore are the product of free 
particle states corresponding to the motion along the xy 
plane and states coming from the confined motion along 
the z axis, which have been calculated by numerical di- 
agonalization of the ID Hamiltonian in the free particle 
basis. We have reached convergence by using the first 
1024 plane waves. We have chosen to investigate in de- 
tail only the pore width oi H = 5.7 A, because this is the 
narrowest slit still presenting a bound state for adsorbed 
hydrogen, according to our model. 

We report in Tab. lIIII the average kinetic and potential 
energies of single hydrogen isotopes confined within the 
H = 5.7 A slit pore. We notice that the average potential 
energies do not change very much with increasing tem- 
perature, a signature of the fact the the first excited state 
is very well separated from the ground state for all the 
three isotopes. The smallest separation is observed for T2 
and is of the order of AiJ/fc^ — 412 K. As a consequence 
the hydrogen isotopes have a very low probability of be- 
ing in the first excited state at the temperatures used in 
this study, and the temperature dependence of the kinetic 
energies reported in Tab. IIIII comes from the different av- 
erage kinetic energies corresponding to the motion along 
the xy plane. 

The difference between the potential and kinetic ener- 
gies of the various isotopes confined within the slit pore 
already indicates that a strong isotope effect is to be 
expected upon adsorption. In fact the the adsorption 
isotherms reported in Fig. [4] show that this is indeed the 
case. 

The isotherms have been calculated using the 
PIGCMC technique and using, for all the isotopes, the 
same number of beads P as in the simulation of adsorp- 
tion within carbon nanotubes. For each pressure we have 
performed 2.5 x 10^ PIGCMC steps for equilibration, fol- 
lowed by 5 X 10^ steps for production. The lateral dimen- 
sions of the pore were such that in saturation condition, 
at least 300 particles were present in the simulation box, 
resulting in a size of about L = 70 A. 
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FIG. 4: Isotherms of adsorption of hydrogen isotopes within 
a slit pore with ii" = 5.7 A at T = 20 K (a) and T = 77 K 
(b). Circles: H2, triangles: D2, squares: T2. The adsorbed 
density is obtained by dividing the total number of molecules 
in the simulation cell by the product of the graphite planes' 
surface and the distance H between them. The lines are a fit 
with a Langmuir type isotherm. 



The computed isotherms for the three isotopes at 
T = 20 K have a similar shape, though they are separated 
by two or three orders of magnitude in pressure; the heav- 
iest specie is adsorbed at lower pressures than the lightest 
one. Moreover the saturation density is markedly differ- 
ent for the three isotopes, with T2 having the largest and 
H2 the smallest. The isotope effect is also sizeable for the 
adsorption at T — 77 K, though not as significant as it 
is at the lowest temperature. The adsorption isotherms 
plotted in Fig. [5] are in good quantitative agreement with 
analogous computer simulation results already appeared 
in the literature . 

Calculation of the selectivity confirms the presence of 
a significant isotope effect. The zero pressure selectivi- 
ties, shown in Table llVi are found to be quite high at 
low temperature, reaching a value of around 1200 for the 
case of T2/H2 mixtures and being of the order of 150 for 
D2/H2 mixtures. The selectivity has a very strong depen- 
dence on the temperature, and drops to values around 4 
and 3 for T2/H2 and D2/H2 mixtures, respectively, at a 
temperature of T = 77 

It is interesting to notice that the values of the zero 
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Adsorbate 


Temperature (K) 


fi/kB (K) 


Kinetic Energy/fes (K) 


Potential Energy/fcs (K) 


H2 


20 


-173.2 


205.4 


-361.8 


D2 


20 


-273.5 


147.7 


-411.3 


T2 


20 


-316.1 


123.0 


-433.2 


H2 


77 


-112.6 


262.4 


-361.7 


D2 


77 


-193.2 


205.1 


-410.9 


T2 


77 


-222.4 


181.1 


-432.2 



TABLE III: Single particle properties of hydrogen isotopes confined in a slit pore with H = 5.7 A, obtained by direct diago- 
nalization of the single particle Hamiltonian. The quantity /i is defined in Eq. (|22|l . 



Mixture 


Temperature (K) 


Zero pressure selectivity 


T2/H2 


20 


1263 


T2/H2 


77 


4.3 


D2/H2 


20 


151 


D2/H2 


77 


2.8 



TABLE IV: Zero pressure selectivities of hydrogen mixtures 
in the H = 5.7 A slit pore 



pressure selectivity calculated for the narrow slit pores, 
although lower than the ones of the narrowest carbon 
nanotubes, are nonetheless larger than those one could 
expect on the basis of geometrical considerations, given 
the 2D confining nature of the carbon nanotubes and the 
ID confining nature of slit pores. In fact, if one assumes 
that the two confining directions are independent the sin- 
gle particle partition function can be written as the prod- 
uct of the partition functions corresponding to the two 
degrees of freedom. As a consequence of Eq. (|24p , then, 
the selectivity could be written as the product of two 
contributions, one for each confining direction. There- 
fore the selectivity of the nanotube should be equal to 
the square of the selectivity of the slit pore, if the two 
systems have the same confinement property along each 
confining direction. 

We have calculated the selectivity as a function of pres- 
sure for adsorption of hydrogen isotopes in slit pores us- 
ing the Boltzmann bias method. We show in Fig. ([5]) and 
([6]) the results obtained at the temperature of T = 20 K 
and r = 77 K respectively. The calculations at each pres- 
sure have been performed by simulating the adsorption 
of mixtures using the Boltzmann bias path integral grand 
canonical method and calculating the finite pressure se- 
lectivity directly from its definition, Eq. (^5)1 . For each 
state point a number between 1.5 x 10^ and 5 x lO'^ of 
Monte Carlo steps were performed for the equilibration 
phase (depending on the pressure) and twice as many 
steps for the production run. In order to have an al- 
most equal number of particles for both species within 
the simulation cell, we set the bulk mole fractions for the 
isotopes to the values indicated in the figures' caption. 

The results for T = 20 K are qualitatively sim- 
ilar to what has already been observed for strongly 
confining nanotubesJ^iii the selectivity rises from its 




Fugacity [bar] 

FIG. 5: Selectivity for hydrogen isotope mixtures as a func- 
tion of pressure in the H = 5.7 A slit pore at T = 20 K. 
The upper panel shows the selectivity of a T2/H2 mixture 
(j/T2 = 8 X 10~*), and the lower panel the selectivity of a 
D2/H2 mixture {y-D2 = 7x 10^'^). The horizontal dashed lines 
correspond to the zero pressure values, reported in Tab. IIVI 



zero pressure value for both tritium/hydrogen and deu- 
terium/hydrogen mixtures, almost doubling its zero pres- 
sure value. This is consistent with the mechanism out- 
lined above in the case of carbon nanotubes. Upon reach- 
ing saturation each adsorbed particle is confined also in 
the xy plane and one would expect the selectivity to raise 
by a factor corresponding to the square of what is ob- 
served in the case of nanotubes. Since the selectivity 
increases by a factor of 1.4 in the case of nanotubes, one 
expects the selectivity in the slit pores to increase by a 
factor of 1.4 x 1.4 ~ 2, as is indeed observed from the 
simulation result at T = 20 K. 

At the higher temperature of T = 77 K the pres- 
sure dependence of the selectivity is less dramatic. One 
can see from Fig. |6] that the selectivity remains almost 
constant over four decades in pressure, possibly show- 
ing a slight tendency to decrease from its zero pres- 
sure value, as is particularly evident in the case of the 
deuterium/hydrogen mixture. In this case the further 
confinement along the xy plane is less efficient than at 
low temperature, and the selectivity remains almost con- 
stant. 
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FIG. 6: Selectivity for hydrogen isotope mixtures as a func- 
tion of pressure in the H — 5.7 A sht pore at T = 77 K. The 
circles shows the selectivity of a T2/H2 mixture (1/T2 = 0.23), 
and triangles the selectivity of a T2/H2 mixture (1/D2 ~ 0.35) 
as a function of the external pressure. The horizontal dashed 
lines correspond to the zero pressure values, reported in 
Tab.HVl 



dependent particles already adsorbed. 

The efficiency of the method comes from the fact that 
one uses as insertion candidates ring polymers already 
equilibrated according to the confining solid-fluid poten- 
tial. We have observed an increase of the acceptance ra- 
tio for insertion moves up to 2 order of magnitudes when 
compared with the standard practice of trying to insert 
ideal gas ring polymer configurations, and a correspond- 
ing decrease of the time needed to perform a simulation. 

The algorithm has been validated by calculating ad- 
sorption isotherms of pure fluids within narrow carbon 
nanotubes and graphite slit pores and comparing them 
with results already published in the literature, or ob- 
tained by us using the standard algorithm. Moreover we 
have calculated the pressure dependence of the isotopic 
selectivity in narrow nanotubes and slit pores and showed 
that even in the case of very strong confinement its value 
increases from the zero pressure limit. 

Although this method becomes progressively inefficient 
near saturation conditions, its simplicity, effectiveness 
and ease of programming makes it a suitable candidate 
for simulating adsorption of quantum gases under condi- 
tions of strong confinement in a wide range of loadings. 

Finally we would like to point out that this simula- 
tion technique can be applied to classical simulation of 
adsorption of small molecules or polymers within narrow 
channels 1^ 



V. CONCLUSIONS 

We have developed an efficient method for the inser- 
tion/deletion moves of ring polymers in path integral 
grand canonical Monte Carlo simulations of fluids ad- 
sorbed in strongly confining geometries. In this case the 
kinetic energy of the adsorbed phase is larger than the 
thermal kinetic energy of a perfect gas at the same tem- 
perature and the standard procedure of inserting ring 
polymers configuration drawn from the ideal gas distri- 
bution is highly inefficient because the polymers are on 
average too large to fit within the narrow channels, re- 
sulting in a very low acceptance ratios. The new algo- 
rithm presented here is based on the idea of taking ring 
polymer configurations from a canonical simulation of in- 

I 

APPENDIX A: PROBABILITY DISTRIBUTION FOR RING POLYMER CONFIGURATION 
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The partition function of a free particle is 

Z = J d^xi {xi I eM-PpV^rn)\xi) = ^ (Al) 
that becomes, after a Trotter expansion, 

Z = / d^xid^X2...(fxp {xi\exp{-(3p^ /2mP)\x2) ■ ■ . {xp\eTq>{-Pp^ /2mP)\xi) (A2) 



p-i 

The integral can be rewritten using the variables defined in Eq. (fT9|) together with the definition Arp = — Ar 

i=l 

Xl - Xp, 



Z = 



p3/2 
"A3~ 



d^rid^Ari . . . d^Arp^i exp 



^ 2 A 

-(Arl + Ar 



Ar^ 



p-i 



Arj.) 
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where K = 27rP/A^. Since the integrand does not depend on ri on can integrate it obtaining 



1 - A'^ 



p3/2 



d^Ari . . . (TArp^i exp 



K 



{Arf + Ar 



Arl_, 



Arl) 



(A4) 



which can be interpreted as the normaUzation condition on the probabihty of observing a ring polymer of P beads 
with separation from the neighboring beads equal to Ari . . . Arp_i, the last one being fixed by the condition of having 
a closed ring polymer. Hence the probability for such a configuration is 



Fri„g(Ari, . . . , Arp_i] 



A^ 



p3/i 



exp 



K 



'Arl + 



Ar^ 



p-i 



Ar?, 



(A5) 
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